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Abstract 


There is a broad consensus that accretion onto supermassive black holes 
and consequent jet formation power the observed emission from ac- 
tive galactic nuclei (AGNs). However, there has been less agreement 
about how jets form in accretion flows, their possible relationship to 
black hole spin, and how they interact with the surrounding medium. 
There have also been theoretical concerns about instabilities in stan- 
dard accretion disk models and lingering discrepancies with observa- 
tional constraints. Despite seemingly successful applications to X-ray 
binaries, the standard accretion disk model faces a growing list of ob- 
servational constraints that challenge its application to AGNs. Theo- 
retical exploration of these questions has become increasingly reliant 
on numerical simulations owing to the dynamic nature of these flows 
and the complex interplay between hydrodynamics, magnetic fields, 
radiation transfer, and curved spacetime. We conclude the following: 


e The advent of general relativistic magnetohydrodynamics 
(MHD) simulations has greatly improved our understanding of 
jet production and its dependence on black hole spin. 


e Simulation results show both disks and jets are sensitive to the 
magnetic flux threading the accretion flow as well as possible 
misaglingment between the angular momentum of the accretion 
flow and the black hole spin. 


e Radiation MHD simulations are providing new insights into 
the stability of luminous accretion flows and highlighting the 
potential importance of radiation viscosity, UV opacity from 
atoms, and spiral density waves in AGNs. 
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Quasar: a luminous 
active galactic 
nucleus that 
outshines its host 
galaxy 


X-ray Binary (XRB): 
X-ray emitting 
source powered by 
accretion from 
companion star onto 
a black hole or 
neutron star 


Active galactic 
nucleus (AGN): 
compact emitting 
region in the center 
of a galaxy thought 
to be powered by 
accretion onto a 
supermassive black 
hole 
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1. INTRODUCTION 


The early observations of quasars and the realization that they must be distant objects 
(Schmidt 1963) already indicated a need for objects much smaller than galaxies to be 
producing immense energies to power their optical emission and radio lobes. Accretion 
of interstellar gas onto a supermassive black hole was identified as the only viable sources 
for such energy (Salpeter 1964; Lynden-Bell 1969). 

The first detailed model of the vertical and radial structure of disks with angular mo- 
mentum transport was provided by Shakura & Sunyaev (1973) with additional insights 
from Lynden-Bell & Pringle (1974) and general relativistic generalization by Novikov & 
Thorne (1973). The key assumptions in the model were that the vertical scale height of the 
disk was small compared with the radius everywhere and that the stress leading to angular 
moment transport could be approximated as a constant fraction (a) of the total pressure in 
the disk midplane. Even 45 years later, this seminal work remains the primary basis of our 
understanding of accretion onto black holes in XRBs and AGN and we refer to this as the 
standard model of accretion disks. Although these steady disk models suffer from instabil- 
ities discussed below, they were successful in broadly explaining the spectral properties of 
XRBs (see e.g. Remillard & McClintock 2006; Done, Gierliriski & Kubota 2007) and the 
optical/UV continuum of luminous AGNs (Shields 1978; Kolykhalov & Sunyaev 1984). 

Although this standard model is thought to apply in the luminous accretion regime, 
there are observational and theoretical arguments that this model breaks down when the 
accretion rate is very low or high when compared with the Eddington (or critical accretion) 
rate, 


_ 4nGMm; 
H NOTC 


Mr (1) 
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where M is the mass of the black hole, 7 is the radiative efficiency, and other values are 
fundamental constants with standard meanings. This is the accretion rate that would 
produce an Eddington luminosity of emission (Le œ 10°°M/Mo) if the radiative efficiency 
remains constant (typically assumed to be around 10%). At low accretion rates (M X 
0.01Me), flows are thought to be hotter and less radiatively efficient (Yuan & Narayan 2014). 
In these models, number densities are low enough that Coulomb collisions are less efficient 
at coupling electrons to ions and radiative cooling is weak, leading to flows dominated by 
radial advection of thermal energy (Narayan & Yi 1994) and/or outflows (Blandford & 
Begelman 1999). At high accretion rates (M > Mg), the assumption of a thin disk is no 
longer self-consistent as the vertical scale height of the inner disk becomes comparable with 
radius. The flow is thought to again become less radiatively efficient as the photon diffusion 
time becomes longer than the radial inflow time, allowing a larger fraction of the dissipated 
energy to be advected into the black hole. Models including these effects are called slim 
disks (Abramowicz et al. 1988). 

All these accretion disk models rely on an anomalous viscosity or stress (o prescription) 
to explain angular momentum transport because realistic physical viscosities are too low 
by many orders of magnitude. Much of the theoretical work in the field for several decades 
focused on discerning the origin of this anomalous stress. After the identification of the 
role of the (MRI; Balbus & Hawley 1991, 1998) in driving magnetohydrodynamic (MHD) 
turbulence in disks, it is now widely believed that MRI turbulent stresses provide the dom- 
inant angular momentum transport mechanism. To date, this has not led to a significant 
revision of the standard model because the Maxwell and Reynolds stresses provided by the 
MRI turbulence are broadly consistent with the a prescription ansatz (Balbus & Papaloizou 
1999). However, in addition to disk internal stresses, disk outflows can also remove the an- 
gular momentum from the accretion disk and affect the disk structure. This is particularly 
true in the case of MHD-driven outflows, where nonlocal stresses are dominant (Blandford 
& Payne 1982). This suggests that the processes of accretion and outflows may be tightly 
coupled. 

Observations indicate that accretion systems produce various types of outflows. Of 
Assisted 
with large-scale magnetic fields, jets can derive their power from the rotation of either the 
black hole via the Blandford & Znajek (1977) effect or the rotation of the accretion disk 
via the Blandford & Payne (1982) mechanism. The jets can cover many decades in radius, 


particular interest are collimated outflows, more commonly referred to as jets. 


spanning from the event horizon of the black hole to the outskirts of the galaxies and 
beyond. Modeling jets is challenging owing to their nonlinear nature that gets exacerbated 
by the poorly understood disk-jet connection near their base and turbulent interactions 
with the ambient medium as they propagate away from the center. 

Although work on the revision and development of analytic models continues, the dy- 
namic and three-dimensional nature of the problem of coupled radiation transfer and MHD 
in general relativistic spacetimes has proven to be a challenge for further analytic progress 
and simpler modeling efforts. The critical role played by dynamic magnetic fields in the 
launching of jets and evolution of accretion disks has led to an increased focus and reliance 
on numerical simulation modeling. In this review, we summarize recent developments and 
the current status of numerical simulations of disks and jets. We review observational re- 
sults that challenge the existing models and motivate the incorporation of new or improved 
treatments of physical processes in state-of-the-art simulations. 
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Eddington 
luminosity: The 
luminosity for which 
outward radiation 
pressure on the 
electrons equals the 
inward gravitational 
pull on the ions 


Magnetorotational 
instability (MRI): 
Instability driven by 
magnetic fields 
interacting with 
shear that provides 
a source of angular 
momentum 
transport 


Magnetohydrodynamics 
(MHD): The 

combined magnetic 
and fluid dynamics 

of conducting 

plasmas 


Blandford-Znajek 
(BZ): Mechanism by 
which magnetic 
fields are twisted by 
black hole spin to 
form an outflow 
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2. CURRENT STATUS OF OBSERVATION AND THEORY 
2.1. Accretion Disk Theory in AGNs 


There are numerous presentations of the standard disk model equations, so we do not 
attempt to reproduce them here, but in our opinion the seminal paper by Shakura & 
Sunyaev (1973) remains the best place to start. For a detailed discussion of how this model 
applies to AGNs, we recommend the textbook by Krolik (1999). Here, we briefly review 
the key results relevant to this work. Broadly speaking, the accretion flow divides into two 
regions: an inner region, in which radiation pressure is the primary support against gravity, 
and an outer region, in which gas pressure dominates. The transition radius is a function 
of disk parameters and is given by 


— E UH (a M VOS REO 
udi "An Oil 0.1 108M, RO paya 


(2) 


where rg — GM/c? is the gravitational radius, and a is the famous parameter that relates 
the accretion stress Zo to the total midplane pressure via Tre = OPen. The radial- and 
black-hole-spin-dependent quantities Rr, Rz, and Rr are defined in chapter 7 of Krolik 
(1999) and represent general relativistic correction factors. Here, Pio: is the sum of gas 
pressure P, and radiation pressure P, because the standard model does not consider the 
impact of magnetic pressure. One can derive different scalings of quantities of interest in 
these inner and outer regions. For brevity, we focus on the scalings in the radiation-pressure 
dominated inner regions and assume electron scattering dominates the opacity. We find 
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where r is radius. The effective temperature Te is an estimate of the surface temperature 
and pmia and Tia are the midplane density and temperature, respectively. The ratio of 
h/r, where h is the disk scale height for the radiation dominated flow is 
0.1 L Rrrg 


h ; 
r maga n OlLg Ra r` 


(See the sidebar titled Thermal Instability in the Standard Accretion Disk Model.) 

Although the standard accretion model approximately succeeds at explaining SEDs of 
AGN, there are a number of observational and theoretical issues that have arisen. Shortly 
after the model was formulated, it was realized that the model is subject to inflow (or 
viscous) and thermal instabilities (Lightman & Eardley 1974; Shakura & Sunyaev 1976 
when radiation pressure exceeds gas pressure in the disk midplane. Inserting the above 


scalings into standard expression for the radiation and gas pressure yields the scaling, 
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Thermal Instability in the Standard Accretion Disk Model 


The thermal instability of radiation pressure-dominated accretion disks (Shakura & Sunyaev 1976) can be 
understood via a simple argument about how the heating and cooling rates in an accretion flow depend on 
temperature (Piran 1978). In an accretion flow, the cooling rate per unit mass Q. from a disk annulus at 
radius r is simply given by the local radiative flux F, divided by the mass surface density (mass per unit 
area) X so that (accounting for both sides of the disk) 


The rate at which energy is dissipated in a shear flow is related to the product of the accretion stresses Tre 
and the shear rdQ./dr. Integrating over height and dividing by X yields a heating rate per unit mass 


E dQ a bhtreQ 
=g [me (Ja x 


Using the standard model equations in the radiation pressure-dominated regime with electron scattering 
opacity yields bo F, x P,/Y and Tre x P,. Since P, = aT*/3 we have: 


ins JE 
Q+ x Sa Q- o p2 
In deriving the standard model, we assume the disk is in thermal equilibrium with Q- = Q+. Now suppose 


that the temperature is perturbed upward. Because departures from thermal equilibrium happen on the 
thermal timescale while changes in the surface density happen on the much longer viscous timescale we 
assume X is constant. An increase in the temperature causes a large increase in the heating rate due to 
the T? dependence and a somewhat more modest change in the cooling rate due to the T^ dependence. So, 
as the temperature is perturbed upward, the heating rate goes up faster than the cooling rate, leading to 
a thermal runaway. The same argument applies in reverse to downward perturbations in temperature, and 
we conclude that radiation pressure-dominated disks are thermally unstable. 


Hence, for the same Eddington ratio, the radiation pressure in at 10?Mo black hole is 
expected to 100 times more dominated by radiation pressure than a black hole with 10Mo. 
'This extreme dominance of radiation pressure in AGNs could suggest that the thermal 
and inflow instabilities could have a larger impact in AGNs. It is generally thought that 
such instabilities will lead to limit cycle behavior in which the disk oscillates between a gas 
dominated lower branch and an advection dominated upper branch (Nayakshin, Rappaport 
& Melia 2000; Janiuk, Czerny & Siemiginowska 2002). The relevant timescales are the 
dynamical time, thermal time, and inflow (viscous) time given by 


1 T 3/2 
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where Q = (GM/r*)'/? is the Keplerian rotation rate and h/r is the aspect ratio of the 
disk scale height to radius. Time-dependent models yield evolution along stable branches 
occurring on the inflow timescale, whereas evolution between branches occurs on the thermal 
timescale. If the disk is as thin as commonly thought, the inflow timescale is very long 
compared to observing timescales for most systems. 

It was shown that the thermal and inflow instabilities could be suppressed in alternative 
stress prescriptions (Sakimoto & Coroniti 1981) where the stress scaled not with the total 
pressure but gas pressure only. However, as we discuss below, evidence from radiation MHD 
simulations of accretion flows suggests the stress does generally scale with the total pressure 
and, thus, provide some evidence that thermal instability is present (Jiang, Stone & Davis 
2013; Mishra et al. 2016). 


2.2. Observational Status: AGNs Versus XRBs 


On the observational side, there are also a number of potential discrepancies. Koratkar & 
Blaes (1999) provide a somewhat dated but still useful review of several of the problems, 
many of which remain unresolved. Some of the issues that have received the most attention 
are the apparent discrepancies between the standard model and sizes of emission regions 
determined by microlensing (e.g Morgan et al. 2010; Blackburne et al. 2011; Morgan et al. 
2018) and continuum reverberation mapping (see e.g. Edelson et al. 2015; Jiang et al. 2017; 
Homayouni et al. 2019). A summary of the microlensing size scale constraints from Morgan 
et al. (2018) is shown in Figure 1. Their best-fit relation suggests that emitting regions are 
about a factor of 3-4 greater than the standard model predicts, which is consistent with the 
time delays inferred from continuum reverberation mapping. 

Variability is also a challenge for disk models. Although the XRB GRS 1915+109 
famously shows variability that looks somewhat like the limit cycle predictions (Nayak- 
shin, Rappaport & Melia 2000), XRBs show no generic evidence of such variability (Done, 
Gierlifiski & Kubota 2007). In AGNs, the predicted viscous timescales are expected to be 
too long (in some cases, thousands of years) to probe observationally. Instead, the observed 
variability poses a different challenge. Surveys and monitoring campaigns have identified a 
large number of “changing look" AGNs, that show large-amplitude variations in the con- 
tinuum or emission lines on timescale of months to years (LaMassa et al. 2015; MacLeod 
et al. 2016), which are difficult to explain with “viscous” evolution timescales (Lawrence 
2018). The question becomes how does one obtain large-amplitude variability of accretion 
disks on such short timescales? 

Others have looked at the spectral variations of the UV continuum of AGNs with mass 
estimates and found them to be poorly fit by spectral models based on the a-disks (see e.g., 
Kolykhalov & Sunyaev 1984; Laor & Netzer 1989; Sincell & Krolik 1997; Bonning et al. 2007; 
Davis, Woo & Blaes 2007). The discrepancy arises in the sense that the model spectra rise 
too rapidly into the far- to extreme-UV, whereas the observed spectra tend to be flat or 
declining in AF, above ~ 1000À. As a typical example, PG 0052+251 is shown in the left 
panel of Figure 2. Here, the comparison relativistic blackbody disk model has been chosen 
to have the best-fit reverberation mapping mass and an accretion rate scaled to match the 
spectrum in the optical band. This has led to the consideration that some other physics is 
at play, with speculations including advection (Netzer & Trakhtenbrot 2014), reprocessing 
close to the inner disk (Lawrence 2012), or outflows (Slone & Netzer 2012; Laor & Davis 
2014). However, this is also the wavelength range in which dust reddening in the AGN host 
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Figure 1: Discrepancies between microlensing measured sizes and predictions from the 
standard model. The solid line gives the best fit relation and the shaded regions shows the 
1c uncertainties. For comparison, the black dot-dashed curve provides the estimate for thin 
accretion disk models with the same mass (assuming L/Le = 1/3, n = 0.1). The observed 
relation is about a factor of 3-4 larger than the standard model prediction. For comparison, 
the solid lines labeled "Last stable orbit" provide the size of the innermost stable circular 
orbit for non-spinning and maximally spinning black holes as a function of mass. Figure 
adapted from Morgan et al. (2018) with permission. 


galaxy might be important, and much of the discrepancy might be accounted for with dust 
reddening (e.g., Capellupo et al. 2015). Unfortunately, the precise form of the reddening 
curve and the amount of extinction remain subjects of debate (Richards et al. 2003; Gaskell 
et al. 2004; Baron et al. 2016). Other significant problems include the lack of evidence for 
Lyman edge features (Shull, Stevans & Danforth 2012), which tends to be prominent (in 
either emission or absorption depending on model parameters) in detailed spectral models 
(Hubeny et al. 2001). 

The overall picture is one in which the standard accretion disk model roughly predicts 
the presence of the optical to UV continuum but seems at odds with a number of obser- 
vational constraints when considered in detail. This situation contrasts with the case of 
XRBs, which the standard model was developed to explain (Shakura & Sunyaev 1973). 
The phenomenology of black hole XRBs can be rather complicated, but there are a number 
of good reviews (see e.g., Remillard & McClintock 2006; Done, Gierlinski & Kubota 2007). 
The spectral evolution is usually described in terms of a number of spectral states, with the 
high/soft (thermal dominant) and low/hard states being the best-characterized states based 
on the differences in their spectral shape, count rates, and variability. The high/soft state 
features a prominent broad thermal continuum component, which is generally assumed to 
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Figure 2: Left panel: Comparison of observed broadband quasar spectral energy distribu- 
tions from PG 00524-251 (Shang et al. 2005) with relativistic multitemperature blackbody 
models chosen to match the reverberation-mapped mass and with accretion rates chosen to 
match the optical emission. The blue X-ray curves show the level and slope change of the 
X-ray emission. Note that there is a turnover in the data in the ultraviolet band around 
3 x 101? Hz (1000A), whereas the model (red dashed curve) continues rising into the extreme 
ultraviolet. This is characteristic of most of the quasars in the PG sample (Davis & Laor 
2011). Right panel: The best-fit model to the black hole X-ray binary source LMC X-3 
(Davis, Done & Blaes 2006). The BeppoSAX data are shown in black, whereas the best-fit 
model is shown as a green curve. The model has two components, with the best-fit absorbed 
relativistic accretion disk model shown in red. The unabsorbed relativistic accretion disk 
model is show in orange. The model provides a good fit to the data with only two free 
parameters in the disk model. 


be emission from a geometrically thin, optically thick accretion disk, whereas the low/hard 
states is dominated by harder X-ray continuum component that is generally thought to be 
produced by inverse Compton scattering of photons in a hot corona above or interior to the 
disk. Although XRBs can have significant variability on shorter timescales, this is mostly 
associated with the coronal component. The high/soft state tends to be characterized by 
relatively low variability, particularly in disk-like component of emission. 

If one focuses solely on the seemingly disk-dominated high/soft state, then the predic- 
tions of the standard model seem to match the data well (Davis et al. 2005; Dunn et al. 
2011). A good example of this, LMC X-3, is shown in the right panel of Figure 2 with 
the best-fit model being a fully relativistic spectrum derived from self-consistent radiation 
transfer calculations through a vertically varying disk structure (Davis et al. 2005). There- 
fore, a key question is why is there this seeming discrepancy in how well standard disk 
models reproduce the properties of AGNs and black hole XRBs? One possibility is that 
we have overstated the agreement between the models and XRB observations. Although 
there are compelling models for the origin of the coronae and the nature of the transitions 
between the different spectral states (e.g., Esin, McClintock & Narayan 1997; Ohsuga et al. 
2009), there is no consensus and this remains an active area of research. Whether or not 
one accepts the premise that standard models work for XRBs and fail for AGNs, there are 
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several reasons to expect that accretion flows in AGNs and XRBs might differ. We focus 
on three possibilities. 

'The first possibility is that the standard model predicts that disks around supermassive 
black holes should be more radiation pressure-dominated than those around stellar mass 
black holes. For accretion rates typically inferred in high/soft state XRBs, these flows 
are expected to be radiation pressure-dominated, but possibly only by factors of ~ 10 
or less. In contrast, the inner radii of the most massive and luminous AGNs could be 
radiation-dominated by factors of a million or more. It is therefore conceivable that the 
impact of thermal and inflow instabilities leave XRBs relatively unchanged but radically 
alter the structure of AGN disks. The main argument against radiation pressure driven 
instabilities reshaping disks is that there is no pervasive evidence of the variability that 
might be expected (see e.g, Done, Gierliriski & Kubota 2007), but it is conceivable that 
such flows find an alternative equilibrium in which such variability is quenched. 

A second possibility is that the much larger opacities present at the lower temperatures 
in AGNs due to dust and ions could significantly alter the disk structure or drive outflows. 
There is considerable evidence of powerful outflows in AGN systems. Broad absorption-line 
quasars show blueshifted absorption lines in the optical band with velocities 5 10,000 km/s 
Weymann et al. 1991) and some sources show ultrafast outflows in their X-ray spectra 
Tombesi et al. 2010). The former are thought to possibly be driven by radiation pressure 
on optical and UV lines (Murray et al. 1995) whereas the latter may also be driven by 
radiation pressure on electrons. Numerical simulations of radiation pressure-driven outflows 
Proga & Kallman 2004) suggest substantial mass loss could be occurring at hundreds of 
gravitational radii or less, where the accretion disk continuum spectrum is formed. Such 
high levels of mass loss could have important implications for the accretion disk spectrum 
Laor & Davis 2014). 

'The simulation of line-driven winds is challenging because of the need to account for 


the impact of X-rays on the ionization state of the gas (Proga & Kallman 2004). Even if 
the X-rays are modeled through radiation transfer, one still needs to use an approximate 
force multiplier scheme to account for the large number of bound-bound transitions that 
may play a role in the acceleration (Castor, Abbott & Klein 1975). In principle, a simpler 
approach uses mean opacities (e.g. Rosseland or Planck mean opacities) to account for the 
atomic species. The Rosseland mean is usually used to specify the radiation force, because 
it is the relevant opacity for a radiation-weighted mean in the diffusion limit (Rybicki & 
Lightman 1986). Opacities in approximately the relevant range of AGN disks have already 
been constructed for computing the evolution of massive stars, and Figure 3 shows the 
Rosseland mean opacities from the OPAL project (Iglesias & Rogers 1996). 

Though it looks like only a modest enhancement, the feature near T ~ 3 x 10° K known 
as the Fe opacity bump is thought to produce extreme variability and help drive outflows 
in the envelopes of massive stars (Jiang et al. 2018). Looking at equations (3) and (4), 
we see that inner regions of near-Eddington accretion rate flows tend to be hotter than 
this in the inner-most radii of the disk, so electron scattering will dominate. However, as 
one moves out, there is always a broad range of radii that are in the relevant temperature 
and density regime. For black hole masses of > 109 Mc this happens in the inner several 
hundred gravitational radii, where the bulk of optical to UV radiation is produced. It is 
conceivable that dynamics driving extreme convection and outflows in massive stars are 
also generically occurring in massive black holes. 

A third set of considerations is the impact of the different environments and feeding of 
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Figure 3: Comparison of a sum of electron scattering opacity with Rosseland mean opacities 
from the OPAL project (black curves) and a Kramer's-like estimate free-free opacity (red 
curves). Line shapes correspond to densities of 107”? (solid), 107° (dotted), and 1078 
(dash-dotted) g/cm?. The modest bump at ~ 3 x 10° K is the Fe opacity bump. 


the two systems. AGNs are thought to be fed predominantly from gas in the large-scale 
interstellar medium (ISM) of their host galaxies, whereas XRBs are fed by Roche lobe 
overflow or wind from a companion star. Outer regions of AGN disks are expected to be 
self-gravitating, and either system could be subject to H ionization instabilities (Lasota 
2001). AGN accretion flows are also expected to be coincident with nuclear star clusters of 
their host galaxies and may have non-trivial interaction with stars or stellar remnants (see, 
e.g., Miralda-Escudé & Kollmeier 2005; McKernan et al. 2012). 

Can these differences in the outer regions of the disks lead to differences in the flow near 
the inner radii, where most of the emission is thought to originate? Because most of the 
key timescales (inflow, thermal, dynamical) decrease with radius, one might expect the disk 
to lose information about its large-scale feeding and slowly adjust to a similar steady state 
model if other considerations (e.g. opacities, radiation pressure) do not drive differences in 
the innermost regions. However, this argument may not hold if the large-scale geometry of 
the magnetic fields is influenced by conditions in the outer flow. Evidence from numerical 
simulations suggests that the presence of a poloidal field may impact both the accretion 
disk and jet formation, but the presence of poloidal field may be partially dependent on 
whether such disks can efficiently accrete a large-scale poloidal field to their inner regions 
(Lubow, Papaloizou & Pringle 1994; Beckwith, Hawley & Krolik 2009; Zhu & Stone 2018) or 
generate it in situ (Beckwith, Hawley & Krolik 2008; McKinney, Tchekhovskoy & Blandford 
2012; Liska, Tchekhovskoy & Quataert 2020). 

We emphasize that radiation plays a key role, particularly in the first two of these con- 
siderations and may in the third as well. If one wants to understand the differences between 
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AGNs and XRBs, then models and simulations in which radiation is treated explicitly are 
required. The physics involved also requires dynamical and multidimensional treatments. 
Such considerations strongly motivate general relativistic radiation MHD numerical simu- 
lation to study the evolution of such flows with realistic opacities, initial conditions, and 
boundary conditions appropriate for AGNs. Of these questions, the most difficult problem 
for simulations to address is the effects of differences in large-scale feeding. The discrep- 
ancies in timescales between the largest and smallest radii are a significant impediment to 
numerical simulations, which must resolve the shortest flow times in the inner regions but 
must be run for many inflow times to achieve steady states at large radii. Nevertheless, 
simulations can examine what impact, if any, the magnetic field strength and geometry 
have on the accretion flow because the presence and strength of the poloidal magnetic field 
can be controlled through the initial conditions. Hence, we think that the key questions 
driving numerical simulation of luminous accretion flows in AGNs are the following: 


e What role does radiation and radiation pressure play in AGNs? Are radiation 
pressure-driven instabilities present, and how do they manifest in the flow structure 
and dynamics? 

e Are a significant fraction of observed AGNs accreting near enough to (or above) the 
Eddington limit so that effects of advection and outflows become important? How 
do these effects modify the light curves and SEDs of accretion flows in this limit? 

e What role do dust and atomic opacities play in the dynamics and structure of AGNs? 
Are there significant outflows of gravitationally bound or unbound material launched 
from the inner regions of AGNs even when the accretion rates are below the Eddington 
limit for electron scattering opacity? 

e What role do magnetic fields play in AGNs? Does magnetic pressure support change 
the structure of the disks? Is there any reason to believe that magnetic field geometries 
differ between AGNs and XRBs? 


Many of these effects are equally important for understanding accretion in XRBs and 
possibly ultraluminous X-ray sources, which may be the result of accretion above the Ed- 
dington limit. However, here we focus on the issues specific to AGNs. 


2.3. Status of the Standard Jet Model Interpretation of AGNs 


It is generally agreed that for a black hole to launch jets it needs to be accreting. However, 
the ingredients necessary for jet formation are poorly understood. Is it the black hole spin 
or the rotation of the accretion disk that powers the jets? Complicating the situation, the 
ability of black hole systems to produce jets depends on the state of the accretion flow: 
'The same black hole can exist in jetted and jet-phobic states depending on how it is fed. 
Perhaps the best evidence for this comes from microquasars: XRBs show radio emission 
that typically varies during spectral state transitions (Fender, Belloni & Gallo 2004). Before 
the transition, the system is in a low/hard accretion state characterized by a nonthermal 
spectrum and accompanied by weak continuous jets. As the transition starts, the luminosity 
of the system goes up, and powerful transient jets emerge. As the luminosity declines, the 
jets disappear, revealing a near-thermal accretion flow. The power of transient jets appears 
to correlate with the black hole spin (Narayan & McClintock 2012; Steiner, McClintock & 
Narayan 2013; but see Russell, Gallo & Fender 2013), suggesting that black hole rotation 
plays a role in powering these outflows. However, how this correlation emerges is unclear. 
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Figure 4: Illustration of jet formation by magnetic fields. (a) Consider a purely poloidal 


De, toroidal field vanishes, B, = 0) field line attached on one end to a stationary “ceiling” 
(which represents the ambient medium and is shown by a hashed horizontal line) and on 
the other end to a perfectly conducting sphere (which represents the central black hole 
or neutron star and is shown by a gray filled circle) rotating at an angular frequency Q. 
(b) After N rotations, at time t — t1, the initially purely poloidal field line develops N 
toroidal loops. This magnetic spring pushes against the ceiling with an effective pressure 
Dm ^v B?, /8r due to the toroidal field, B,. As time goes on, more toroidal loops form, and 
the toroidal field becomes stronger. (c) At some later time, t = t2, the pressure becomes so 
large that the magnetic spring, which was twisted by the rotation of the sphere, pushes away 
the ceiling and accelerates the plasma attached to it along the rotation axis, forming a jet. 
Asymptotically far from the center, the toroidal field is the dominant field component and 
determines the dynamics of the jet. (d) It is convenient to think of the jet as a collection 
of toroidal field loops that slide down the poloidal field lines and accelerate along the jet 
under the action of their own pressure gradient and tension (hoop stress). The rotation of 
the sphere continuously twists the poloidal field into new toroidal loops at a rate that, in 
steady state, balances the rate at which the loops move downstream. The power of the jet is 
determined by two parameters (equation 11): the rotational frequency of the central object, 
Q, and the radial magnetic flux threading the object, @. Figure adapted from Tchekhovskoy 
(2015). 


To understand the disk-jet connection, let us first review the basics of jet formation 
without attempting to connect it to the disk. For simplicity, let us consider a perfectly 
conducting spinning sphere, shown in Fig. 4(a), which is meant to represent the central 
compact object (a black hole, neutron star, white dwarf, or even normal star), and a per- 
fectly conducting “ceiling”, which is meant to represent the ambient medium. Suppose a 
magnetic field line connects the sphere to the ceiling. As the sphere rotates, it coils up the 
field line into a magnetic spring (Fig. 4b), which pushes the ceiling away and accelerates 
under the action of its own pressure (Fig. 4c). One can view this acceleration in an alterna- 
tive way (Fig. 4d): the rotation of the central sphere continuously reprocesses the initially 
vertical field line into toroidal field loops that emanate from the sphere and slide up the 
jet. As they do so, they expand, their pressure drops, and the pressure gradient accelerates 
them away. The situation for black holes, qualitatively, is rather similar. Even though a 
black hole does not have a physical surface, black hole rotation leads to the dragging of 
the inertial frames near the black hole. This causes the magnetic field lines to rotate in a 
surprisingly similar fashion to the perfectly conducting sphere. 
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Blandford & Znajek (1977) showed, in the limit of slow rotation, that a spinning black 
hole immersed into a large-scale vertical magnetic field would produce jets of power 


Paz = k(a/2r,)@Zuc, (11) 


where ®puy is the magnetic flux threading the black hole event horizon, —1 <a < lisa 
dimensionless black hole spin parameter, and k is a dimensionless proportionality factor. 
Ignoring constant prefactors, we can obtain this expression from dimensional analysis by 
writing that the jet power is the product of magnetic energy density, x B?, the cross-section 
of the base of the jet, x ta and the speed v ~ c with which the energy flows through the 
jets, giving P « a?r? Ba. Here, we introduce the a? prefactor to account for the variation 
of jet power with black hole spin (it has to be an even power of spin because spin sign 
change, by symmetry, leaves the power unchanged). Switching from the field strength, B, 
to the magnetic flux, Ọga ~ Br, gives P œ (a/rg)’?®yc, which has the same scaling as 
in eq. (11). 

With the advent of numerical simulations, it became clear that though equation (11) 
works well for small values of a, it underestimates the jet power for rapidly spinning black 
holes (Komissarov 2001). Comparison against numerical solutions shows that replacing ry 
with the event horizon radius, rg = rg[1 + (1 — a?)'/?], gives a more accurate expression, 


Pret = k(a/ru) ®pue x f(a), (12) 


that serves well for most practical purposes at a < 0.95 with f(a) = 1 (Tchekhovskoy, 
Narayan & McKinney 2010). A higher-order correction, f(a) = 1 + 0.35(arg/rg)? — 
0.58(ar,g/ri)*, allows eq. (12) to maintain accuracy all the way up to a = 1 (Tchekhovskoy, 
Narayan & McKinney 2010; Pan & Yu 2015). 

Thus, there appears to be a clear relationship between the black hole spin and jet power. 
What, then, causes the jet power to change for any given black hole? Since in AGNs the 
black hole spin does not change in our lifetime, changes in jet power can only come from the 
variations in the black hole magnetic flux, fgg. To understand the disk-jet connection, we 
therefore need to understand how the accretion physics determines the value of ®gxy. If we 
were to remove the accretion disk, then the black hole would lose its magnetic flux: By the 
no hair theorem, the black hole can only have three hairs — mass, spin, and charge (Misner, 
Thorne & Wheeler 1973). Thus, the accretion disk holds the magnetic flux on the black 
hole and prevents it from slipping away: Pressure of the magnetic flux on the black hole 
must be in some way balanced by the pressure of the accretion flow. Here, we can again 
turn to the dimensional analysis and write that the magnetic pressure on the black hole is 
proportional to 62,;. This balances the ram pressure in the inflow, which should be roughly 
proportional to the mass accretion rate, M. To characterize the relative strength of the two, 
we introduce a dimensionless black hole magnetic flux, dn = ®pu/ (M r2c)/ ? Its value is 
set by the interaction with the disk. The nonlinearity and intrinsic time-variability of the 
disk-jet interaction makes numerical simulations an attractive approach for constraining 
the allowed range of values of dn. As we discuss in Section 5, the allowed range of dn 
spans the range from zero to a maximum value around 50 for which the black hole magnetic 
flux becomes as strong as the gravity that keeps the disk on an orbit around the black hole 
(Tchekhovskoy, Narayan & McKinney 2011). In this magnetically arrested disk (MAD) 
state, the magnetic flux is as strong as possible and weakly depends on the black hole spin, 
gmap 2 70(1 — 0.38arg/re)h3/2, where h = r x 0.3ho0.3 is the half-thickness of the disk. 
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This corresponds to the following maximum jet power: 
Puan = "Map Mc? © 1.3ho,3a Mc, (13) 


where we used the fact that the jet efficiency is nmap = késjApa (rg /ru)? f (a) = 1.3a? ho. 
(Tchekhovskoy 2015). Thus, the jet power Piet takes on a value between 0 and Dan 
depending on the strength of the black hole magnetic flux py relative to its maximum 
possible value, maD: 


2 
0 € Pet = ( $BH ) Puan < Puman = 1.3h0.307Mc’. (14) 
ỌMAD 
Now that we understand the jet power, there are still many questions left unanswered. 
These questions are driving detailed studies in various branches of astrophysics include the 
following: 


e Where does the magnetic flux powering jets come from? Does it need to be dragged 
inward from the ISM or can it be generated inside the disk in situ? 

e How are transient jets launched during spectral state transitions in XRBs? Does an 
analog of the state transitions and transient jets exist in AGNs? 

e How do jets accelerate to relativistic velocities? Is radiation pressure important in 
launching and accelerating the jets? 

e Do jets heat up the interstellar gas and affect galaxy evolution? Can jet feedback 
lead to the M - ø relationship? 

e What makes jets shine? What can we learn from the observations of the black hole 
shadow with the Event Horizon Telescope? What can the jets tell us about the 
strong-field gravity and general relativistic frame dragging that birthed them? 


We discuss several of these questions in Sections 3 and 5. 


3. OVERVIEW OF ACCRETION DISK AND JET SIMULATION METHODS 
3.1. General Relativistic MHD 


'The numerical simulation of accretion flows onto black holes has a long history, with Wilson 
(1972) already considering the flow of fluid with nonzero angular momentum in the Kerr 
spacetime. Early work primarily focused on hydrodynamics models (e.g., Hawley, Smarr 
& Wilson 1984), but the realization that the MRI could provide the necessary angular 
momentum transport led to the development and application of MHD methods (Balbus & 
Hawley 1998). With this realization it became standard for numerical simulations to solve 
some version of the MHD equations, which includes conservation of mass, 


Op 
ale : = 1 
ay +V.: (ov) =0, (15) 
conservation of momentum, 
oro) +V-(pvv—- BB + P*) = -pVo, (16) 
conservation of energy, 
OE P 
— +V. |(E + P*)v— B(B:v)| 2 pv. Vo, (17) 


Ot 
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Numerical Methods for Solving the MHD Equations 


There are a significant number of different methods available for simulating MHD flows and an even larger 
number of codes that implement them. Most methods for solving fluid dynamics equations fall into two 
categories: Eulerian or Lagrangian. Eulerian methods tend to solve finite difference representations of the 
equations on a fixed mesh. Finite volume (Godunov) methods (LeVeque 2002) are particularly popular 
for their shock-capturing capabilities because shocks commonly arise in astrophysical flows. Lagrangian 
methods tend to utilize particle-based representations with fluid properties advected along with the particles, 
such as smoothed particle hydrodynamics (Liu & Liu 2003). Other examples include moving mesh codes that 
carry a time variable mesh along with the particles or fluid flow (e.g., Springel 2010; Duffell & MacFadyen 
2011). 

Both Eulerian and Lagrangian codes have their merits, but the majority of the simulations discussed 
in this review utilize mesh-based Eulerian schemes. This is partly due to their shock-capturing capabilities 
and their simplicity relative to moving mesh methods. Another major consideration is the need to preserve 
the divergence-free nature of the magnetic field: V - B = 0. Although it should always hold in nature, 
this condition is not preserved by all integration methods. Several codes address this issue by allowing the 
divergence to develop but then limiting it via a divergence cleaning scheme. Other integration schemes are 
specifically designed to preserve the divergence-free condition to machine precision. A popular example is 
the constrained transport scheme (Evans & Hawley 1988), which is most easily implemented on structured 
Eulerian mesh (c.f. Mocz, Vogelsberger & Hernquist 2014). 


and the induction equation, 


OB 


Here, p, B, v are density, magnetic field and flow velocity, P* = (P; + B?/2)1 (with | the 
unit tensor), P, is the gas pressure, and the magnetic permeability u = 1. The total gas 


energy density is 


B= E, + lps (19) 


and P is the gravitational potential. The precise form of these equations varies depending on 


application. Sometimes the energy equation is omitted and an isothermal equation of state 


is adopted. Additional source terms may be added to the right hand sides of the equations 


of momentum and energy, such as the effects of radiative cooling, heating, or radiation 


forces. This may require the solution of additional equations, such as the radiation transfer 


equation or its moments. 


Early MHD simulations of accretion focused on shearing box simulations (Hawley, Gam- 


mie & Balbus 1995; Brandenburg et al. 1995; Stone et al. 1996), which adopt shearing 
boundary conditions and add the Coriolis force to mimic the effects of fluid rotation in 


an accretion disk but only simulate a small patch of the disk. Although these simulations 


are still widely used to study the structure, dynamics and thermodynamics of the disk at 


high resolution, concerns about convergence with resolution (Pessah, Chan & Psaltis 2007; 
Fromang & Papaloizou 2007; Ryan et al. 2017) and dependence on box size (e.g. Simon, 
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Beckwith & Armitage 2012; Shi, Stone & Huang 2016) remain. (See the sidebar titled 
Numerical Methods for Solving the MHD Equations.) 

Global simulations of MHD black hole accretion flows were first performed with pseudo- 
Newtonian potentials (Hawley & Balbus 2002) but these were supplanted by general rel- 
ativistic MHD (GRMHD) simulations of accretion flows onto spinning and nonspinning 
black holes (Gammie, McKinney & Tóth 2003; De Villiers, Hawley & Krolik 2003). In 
GRMHD, the covariant generalizations of the MHD equations are conservation of mass, 
stress-energy, and the relativistic Maxwell’s equations, which are evolved on a choice of 
coordinates (usually Boyer-Lindquist or Kerr-Schild) dictated by the Kerr spacetime. 
The first GRMHD simulations lacked realistic cooling and, therefore, best approximate 
low-luminosity accretion flows. Nevertheless, they have led to a much better understanding 
of flow in the innermost regions near the black hole, including accretion in the plunging 
region, driving of outflows, and launching of jets. 

In particular, it became clear that in the presence of large-scale vertical magnetic flux 
jets are a typical outcome of radiatively inefficient, geometrically thick black hole accretion 
disks with h/r ~ 0.3—1 (McKinney & Gammie 2004; De Villiers et al. 2005; McKinney 2005; 
Hawley & Krolik 2006). Radiatively efficient, geometrically thin disks with h/r < 0.05 were 
simulated via the inclusion of cooling terms that kept the disk thickness at a desired value 
(Shafee et al. 2008; Penna et al. 2010; Noble, Krolik & Hawley 2010). These simulations did 
not produce jets, in agreement with the observations that show no detectable radio emission 
(and, hence jet activity) from thin disks in XRBs and the vast majority (9096) of quasars. 
The remaining 1096 of quasars show powerful radio jets whose origins are still a mystery. 
Recent GRMHD simulations show that even very thin disks, h/r ~ 0.015, are capable of 
producing powerful jets and provide potential resolution to this puzzle (see Sec. 5.1 and 
Liska et al. 2019c). 

Interestingly, it appears that higher resolutions than typical (e.g., 256? cells) are required 
to achieve convergence in the value of the effective a-viscosity parameter in global GRMHD 
simulations of weakly magnetized accretion disks (Porth et al. 2019). Global simulations 
of strongly magnetized accretion disks are less sensitive to resolution (White, Stone & 
Quataert 2019). 


3.2. Radiation Hydrodynamics 


Adding and improving radiation transfer treatments in numerical simulations of accretion 
flows is an important focus of recent work. The relevant equation to solve is the radiative 


transfer equation, 
10L, 
c Ot 


+n- VL =h - wl, (20) 


where I, is the specific intensity, rj, is the emissivity, a, = &,p is the extinction coefficient, 
and ky is the opacity. The form on the right-hand side is general if extinction includes 
scattering opacity contributions and the emissivity accounts for scattered radiation. Because 
l, is a function of position, angle, and frequency, it is usually computationally expensive 
compared to the standard MHD equations, which only depend on position. For this reason, 
some treatments focus on its angle and frequency-integrated moments corresponding to 
conservation of radiation energy, 


OE, 
Ot 


+V-F, = cp (kpaT* — kgEr) 3 (21) 
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and conservation of radiation momentum, 
1 OF, 
c? Ot 


Here, Er, Fr, and P, are the radiation energy density, flux, and pressure tensor, respectively. 


Lob "Cp, (22) 


The opacities kp, Kz, and xp correspond to the Planck, energy, and flux mean opacities. 
In most applications, the Rosseland mean is used for &r and the Planck mean is used for 
KE. 

For the sake of brevity, we will summarize some of the most salient points and refer 
the reader to Mihalas & Mihalas (1984) for a comprehensive discussion of theses equations. 
Although these equations look simple, one should note that, as written, all variables are in 
the Eulerian frame, but the simple isotropic forms for the emissivities and opacities only 
hold in the comoving frame. Hence, one usually must Lorentz transform the source terms, 
expand the right-hand sides in powers of v/c, or use a comoving frame approach which 
introduces additional acceleration terms. Compton scattering also generally introduces 
additional terms. The source terms on the right hand side of equations (21) and (22) 
represent the transfer of energy and momentum from the radiation field to fluid. Hence, 
the negative of these terms are added to the right-hand side of equations (16) and (17), 
corresponding to net heating/cooling by the radiation field and the radiation force. 

Early efforts to include radiation predominantly focused on simulations utilizing flux- 
limited diffusion (FLD; Levermore & Pomraning 1981). In FLD, only equation (21) is 
retained, and F, is computed from the gradient of E. using a diffusion approximation. A 
limiter is utilized to keep |F,| € cE, in the optically thin regime. This method has been used 
in the context of shearing box (Turner 2004) and global accretion flow (Ohsuga et al. 2005) 
simulations. More recently, several groups have developed general relativistic radiation 
transfer modules based on two moment methods with M1 closure (Sadowski et al. 2013; 
McKinney et al. 2014; Mishra et al. 2016). 
or their relativistic generalization (corresponding to conservation of radiation stress-energy) 


In M1 closure schemes, equations (21) and (22) 


are solved. Because P, appears in equation (22), it must be specified. This is generally done 
by introducing a new tensor called the Eddington tensor f = P,/ E, which characterizes the 
angular distribution of the radiation field. In the M1 scheme, f is computed as a function 
of E, and F,. 

Other approaches include direct solution of the angle-dependent transfer equation 
(Jiang, Stone & Davis 2014b) or the variable Eddington tensor (VET) method (Jiang, 
Stone & Davis 2012). In the direct solution approach, a frequency-averaged version of 
equation (20) is solved explicitly for a fixed number of angles in each zone. In the VET 
method, equations (21) and (22) are integrated directly, but a time-independent solution 
of the transfer equations is used to compute the Eddington tensor. Becasue both the FLD 
and M1 methods make assumptions about the angular distribution of the radiation field by 
prescribing relations among Er, F, or P,, methods involving direct solution of the transfer 
equation are generally thought to be more accurate. The trade-off is that they are also more 
computationally expensive and are algorithmically more complex to implement. For this 
reason, incorporating general relativity in the radiation transfer equation is significantly 
more challenging, so simulations performed with these methods have all been nonrelativis- 
tic MHD-based, whereas many of the simulations performed with M1 methods have been 
general relativistic MHD-based. The formalism for the general relativistic transfer equation 
is well-developed (Cardall, Endeve & Mezzacappa 2013) and efforts to implement it in black 
hole accretion simulations are ongoing. 
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4. MHD SIMULATIONS OF AGN ACCRETION DISKS 


'To date, there have been relatively few simulations directly aimed at studying supermassive 
black hole accretion disks. Prior to the addition of radiation transfer in simulations, there 
was usually no explicit dependence on the black hole mass in the simulations. This is 
because the bare MHD equations can be arbitrarily rescaled by an explicit choice of length 
scale or characteristic density. Choosing the characteristic length scale specified the mass, 
and choosing the density then specified the mass accretion rate. Hence, one could scale 
the simulation to the supermassive black hole regime, but there was nothing to distinguish 
the dynamics or thermodynamics of the rescaled simulations from that of stellar mass 
black holes. However, this rescaling freedom is lost when radiation is added, because the 
dependence of opacities and emissivities on temperature and density enforces an explicit 
choice of length and density scales. In this sense, radiation transfer is fundamental to the 
distinction between supermassive and stellar mass black holes. Although we have learned 
a great deal about accretion flows from pure MHD simulations, they cannot tell us about 
the differences between AGN flows and XRB flows, so we choose to focus on radiative 
simulations in this section of the review. 

The first sets of radiation MHD simulations of accretion flows utilized the FLD approxi- 
mation (Levermore & Pomraning 1981; Turner & Stone 2001) to perform local shearing box 
simulations (e.g. Turner et al. 2003; Turner 2004; Hirose, Krolik & Stone 2006). With the 
exception of Turner (2004), these simulations focused on patches of disk around ~ 10M, 
black holes, where the radiation-to-gas pressure ratio tended to be lower. At the same time, 
efforts were underway to use FLD to study the global structure of accretion flows. Initially, 
this took the form of two-dimensional, viscous hydrodynamic flows (Ohsuga et al. 2005), 
but was eventually generalized to MHD treatments (Ohsuga & Mineshige 2011). These 
calculations were mainly aimed at studying the variation of the flow with accretion rate 
and also focused on 10Mọ black holes. 


4.1. Thermal Instability in Radiation Pressure Dominated Accretion 


These radiative MHD shearing box simulations were particularly interesting for testing the 
predictions of thermal instability in radiation pressure-dominated accretion flows (Shakura 
& Sunyaev 1976). The initial results from the Zeus FLD simulations confirmed that gas 
pressure-dominated disks were stable (Hirose, Krolik & Stone 2006), as expected. More 
intriguingly, simulations found evidence for thermal stability even when radiation pressure 
dominates gas pressure (Turner 2004; Hirose, Krolik & Blaes 2009). However, these results 
were contradicted by Athena simulations utilizing the VET method (Jiang, Stone & Davis 
2013), which generally found runaway heating or collapse. Simulations using the Athena 
FLD module also showed indications of thermal runaway, so the discrepancy was not purely 
a result of the radiation transfer algorithm. Instead, Jiang, Stone & Davis (2013) attribute 
the discrepancy to the radial width of the Zeus simulation box being too small, finding 
that the simulations became more stable as they narrowed the radial dimensions of the 
simulation domain. 

Despite reaching different conclusions about stability, these simulations all found that 
the disks could survive for many thermal timescales, which is indicative of a possible robust- 
ness against the standard thermal instability. Although the relation of the time-averaged 
stress in MHD simulations is roughly consistent with the a-disk assumption, the shearing 
boxes differed from the standard model assumptions in important ways. Magnetic pressure 
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is not included in standard disk models, but in simulations it plays an important role in 
supporting the surface regions against gravity. Thermodynamics is more complicated in 
MHD runs, with a greater fraction of dissipation occurring near or above the photosphere, 
leading the density to drop faster with height than in standard disk models. And the time 
dependence of the relation between central pressure and stress is less direct, so that in- 
creases (decreases) in radiation pressure did not immediately result in increases (decreases) 
in the scale height or dissipation (Hirose, Krolik & Blaes 2009). 

Turner (2004) was the sole simulation among those above that was run for supermassive 
black hole conditions, but it only included free-free (bremsstrahlung) and electron-scattering 
opacity. Jiang, Davis & Stone (2016) explored the role of opacity in such environments, 
using the Athena VET module to simulate several shearing boxes for conditions around a 5x 
105 M black hole while including OPAL opacities. Despite much higher ratios of radiation 
to gas pressure compared with the lower-mass black holes discussed by Jiang, Stone & Davis 
(2013), the resulting simulations appeared to be thermally stable. In contrast, simulations 
performed with the same parameters, but only including free-free and electron-scattering 
opacity, all showed runaway collapse on very short timescales, much faster than in the lower 
black hole mass runs. Jiang, Davis & Stone (2016) attribute the enhancement in stability 
to a combination of two effects. First, UV opacities introduce an anticorrelation in the 
optical depth to the midplane and central temperature, opposite to the electron scattering- 
dominated assumption. Second, the opacity drives convection that isn't present in the 
previous simulations, which enhances the nonradiative fluxes carried in the simulations. 
However, the dependence of shearing box simulations on box size, combined with time- 
dependent global modeling (Grzedzielski, Janiuk & Czerny 2017) strongly motivates further 
study of UV opacity effects in global numerical simulations. 


4.2. Global Simulations of AGN Accretions Flows 


Local (shearing box) simulations are useful for studying accretion physics at high resolution, 
but one ultimately requires global simulations to obtain a more complete picture of the 
accretion flow. Furthermore, the apparent dependence of thermal instability on radial width 
(or aspect ratio) in shearing box simulations motivates global simulations for which this 
artificial constraint is removed. As already noted, global simulations of black hole accretion 
have advanced significantly over the past two decades, advancing from pseudo-Newtonian 
MHD (Hawley, Gammie & Balbus 1995) to full GRMHD (Gammie, McKinney & Tóth 
2003), sometimes employing ad hoc cooling function to maintain a thin disk (e.g., Penna 
et al. 2010; Noble, Krolik & Hawley 2010). Without radiation transfer and realistic opacities, 
such simulations can be scaled to arbitrary masses. Therefore, these simulations apply to 
both the solar mass and the supermassive black hole regime. However, in practice, such 
simulations have been primarily compared with XRBs (e.g. Kulkarni et al. 2011; Schnittman, 
Krolik & Noble 2013). 

For these reasons, AGN specific MHD numerical simulations are relatively recent and 
relatively few in number. Nevertheless, there are several lines of related inquiry that are 
worth summarizing. Although their simulation was calibrated for solar mass black holes, 
Mishra et al. (2016) performed the first radiation GRMHD simulation of radiation pressure- 
dominated global disk. Consistent with the earlier shearing box simulations results, they 
found a runaway collapse of the disk, approximately on the thermal timescale. In contrast, 
Sadowski (2016) performed global radiation GRMHD simulations in a similar regime that 
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was stable to thermal and inflow instabilities. 

Thus, we have two sub-Eddington accretion rate simulations that provide different pre- 
dictions about the stability of accretion flows in this radiation pressure dominated regime. 
Sadowski (2016) attribute the stability of their simulation to magnetic pressure from toroidal 
magnetic fields, which is broadly consistent with predictions of earlier analytic models that 
suggested magnetic support could impact stability (Begelman & Pringle 2007; Begelman & 
Silk 2017). The difference between the Sadowski (2016) and Mishra et al. (2016) results can 
then be attributed to the magnetic field geometry. The required toroidal fields are produced 
by threading the disk with a large poloidal flux that then gets wound up by the shear in 
the accretion flow; this is a picture that is broadly supported by other recent simulation 
results (Mishra et al. 2020), including some in the AGN regime discussed below (Jiang et al. 
2019). These results suggest that magnetic support can stabilize accretion flows, as long 
as sufficiently large poloidal magnetic fields can be accreted into the inner regions of the 
accretion disk. 

Radiation hydrodynamic simulations of line-driven AGN outflows (e.g., Proga, Stone & 
Kallman 2000) have also been used to model properties of broad absorption line quasars. 
Because these calculations treat the accretion disk and its emission as a fixed boundary con- 
dition, they cannot evolve the disk self-consistently, but it is notable that some simulations 
have mass outflow rates that are a substantial fraction of the assumed inflow rate. Taken 
together, these simulations suggest that the impact of atomic opacities and the enhanced 
radiation pressure in AGNs may lead to distinct differences from XRB accretion flows. 

Radiation GRMHD simulations of the radiatively inefficient accretion flow regime have 
also advanced significantly in the past several years, motivated by the imaging of Sgr A* 
and M87 with the Event Horizon Telescope (Ryan et al. 2018; Event Horizon Telescope 
Collaboration et al. 2019). Because such flows are relatively optically thin, they are well 
suited to Monte Carlo-based algorithms for radiation transfer, such as the BHLIGHT code 
(Ryan, Dolence & Gammie 2015). However, the standard implementation of these algo- 
rithms has not traditionally scaled well to the optically thick regime of luminous accretion 
flows, so these techniques have not been applied to luminous AGNs. 

'To date, only a few radiation MHD simulations have been performed with AGN masses 
and UV appropriate opacities. These AGN-specific simulations have been performed for 
super-Eddington (Jiang, Stone & Davis 2019) and sub-Eddington (Jiang et al. 2019) regimes 
using the Athena--- radiation MHD code with a pseudo-Newtonian potential. The simu- 
lations are initialized with a torus centered at 50-80 gravitational radii, seeded with a weak 
magnetic field. The torus is unstable to the MRI, generating turbulence that allows an 
accretion flow to self-consistently form at radii interior to the torus. 

The generic properties of the accretion flow are broadly similar to the two-dimensional 
viscous radiation hydrodynamic simulations (Ohsuga et al. 2005; Kitaki et al. 2018) as well 
as radiation GRMHD simulations performed for solar mass black holes (McKinney, Dai & 
Avara 2015; Sadowski et al. 2015; Sadowski & Narayan 2016) or super-Eddington accretion 
in tidal disruption events (TDEs; Dai et al. 2018). At the highest accretion rates, the 
accretion disk is geometrically thick, with strong radiation pressure-dominated outflows at 
velocities of 10-3096 of the speed of light. However, in contrast to previous simulations and 
slim accretion disk models, the radiative efficiencies can remain relatively high, as much as 
~ 5%. The highest efficiencies are consistent with earlier Athena simulations in cylindrical 
geometry (Jiang, Stone & Davis 2014a) around solar mass black holes, which attributed 
the enhanced radiative efficiency to vertical advection of energy by the MHD turbulence. 
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At higher accretion rates (~ 500Lgaa/c?), radiative efficiencies drop to ~ 1%, in better 
agreement with slim disk models and the radiation GRMHD simulations of Sadowski & 
Narayan (2016). Kinetic efficiencies are generally lower than in those GRMHD simulations, 
but this may be attributed (at least in part) to the use of the pseudo-Newtonian potential 
and the lack of a BZ process driven by a spinning black hole. 

Perhaps the most intriguing results from the Jiang, Stone & Davis (2019) simulations 
is the excitation of spiral density waves and their contribution to driving accretion via 
Reynolds stresses, which in some (but not necessarily all) cases may be the dominant 
contribution. This enhanced role of density waves is attributed to radiative damping of 
MRI turbulence in the simulations, although the details of the excitation mechanism are 
not yet well understood. Radiation viscosity (Masaki 1971; Mihalas & Mihalas 1984) is 
also more important than previous thought. Radiation viscosity acts on shear flows in the 
same way as standard viscosity but with photons providing the momentum exchange across 
the shear layer. In the super-Eddington simulations, it plays a surprisingly significant but 
ultimately subdominant role in the overall transport of angular momentum. The large 
radiation energy density enhances the role of radiation viscosity because more photons are 
present to transfer momentum. However, the large optical depths and correspondingly short 
mean free paths between electron scattering means there are only small changes in fluid 
velocity differences between scatterings. 

The sub-Eddington simulations (Jiang et al. 2019) offer further surprises. Despite having 
accretion rates as low as 796 of the Eddington accretion rate, these simulations appear to 
reach a steady state out to ~ 20 gravitational radii and last for a duration of many thermal 
timescales without any evidence for thermal or inflow instabilities. As shown by Sadowski 
(2016), the stability is attributed to magnetic pressure support. Although magnetic pressure 
is lower than radiation pressure, it has a smaller scale height and provides the dominant 
support of the gas against tidal gravity. The radiation pressure scale height is large because 
radiation pressure within the disk is roughly constant, due to the dissipation primarily 
occurring in surface layers with relatively low optical depth. Such magnetic pressure support 
has been proposed to aid in explaining the size and timescale discrepancies discussed in 
Section 2.1 (Dexter & Begelman 2019). 

'The large dissipation in surface layers is the result of strong radiative viscosity. In the 
simulation with 796 of the Eddington accretion rate, the radiation viscosity is actually the 
dominant mechanism for angular momentum transport, and more mass flows inward in the 
surface layers than in the midplane of the disk where MRI turbulence dominates. The large 
radiation viscosity is the result of a combination of large photon mean free paths in these 
surface regions and large radiation energy density. 

As shown in Figure 5 the large surface dissipation leads to strong temperature inver- 
sions above the disk, reminiscent of models of AGN coronae. The simulations have both 
temperatures high-enough to potentially produce the hard X-ray coronae and more mod- 
erate temperatures needed to explain the presence of soft X-ray excesses observed in many 
AGNs (Walter & Fink 1993; Kubota & Done 2018). Intriguingly, the fraction of emission 
in optically thin regions is higher in the lower accretion rate simulations, leading to higher 
temperatures consistent with the observed relation that optical to X-ray spectral indexes 
get harder as luminosity decreases (Steffen et al. 2006). Although these results are promis- 
ing, we must keep in mind that all these effects are occurring in close vicinity to the black 
hole, where relativistic effects are important, so we must ultimately test these conclusions 
with general relativistic calculations. 
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Figure 5: Time and azimuthally averaged spatial distributions of gas temperature Ty, for 
the two simulations runs onto 5 x 10? black holes from Jiang et al. (2019). The left and 
right panels have accretion rates of 796 (AGN0.07) and 2096 (AGNO.2), respectively. The 
gas temperature is ~ 10? — 2 x 10? K in the optically thick part of the disk but rapidly 
increases to 10° - 10? K in the optically thin coronal regions. Note how the corona is more 
radially extended in the lower accretion rate simulation than in the higher accretion-rate 
case. 


Although OPAL opacities are used in the Athena--4- simulations of super-Eddington 
and sub-Eddington disks, they have relatively little impact on the resulting dynamics. This 
is because the inner regions of the accretion rate, where the flow reaches a steady state, are 
too hot for there to be significant enhancements in the opacity above electron scattering. 
Future simulations need to focus on larger radii to see the effects of such opacities. 


5. SIMULATIONS OF AGN JETS 
5.1. The Disk-jet Connection 


Numerical simulations are a powerful tool for quantifying the power of relativistic jets and 
understanding the disk-jet connection. A particularly useful way of doing so is through 
measuring jet energy efficiency, or the ratio of jet-to-accretion power, 7 = Piet / Mc. Nu- 
merical simulations by different groups found vastly different values of n: For instance, for 
a = 0.99, the efficiency ranged from ~ 3% (McKinney 2005) to ~ 15% (Hawley & Kro- 
lik 2006). Perhaps not surprisingly, the simulations sample the large parameter space of 
allowed jet powers, which can range from zero power (no jet) to some maximum power. 
In the context of magnetically powered jets, this would map into a range from zero large- 
scale magnetic flux to some maximum value that could be measured by the simulations 
(see eq. 14). To determine the maximum jet power, Tchekhovskoy, Narayan & McKinney 
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(2011) initialized GRMHD simulations with what was then considered an unusually large 
gas reservoir threaded with vertical magnetic flux (Figure 6a): The initial gas torus size of 
20, 000rg was much larger than the typical tori of  50r,. The large torus size translated 
into a great amount of vertical magnetic flux contained in the system, thereby flooding the 
black hole with the magnetic flux and determining the maximum jet efficiency achievable 
for a given spin. As the simulation started, the MRI led to the accumulation of gas and 
magnetic flux on the black hole (Figure 6b). This resulted in black hole magnetic flux 
increasing (Figure 6f) to the point that it overcame the force of gravity that keeps the gas 
on an orbit around the black hole: The flux became strong enough to escape from the black 
hole by tearing through the disk (Figure 6c,d). Such periods of flux expulsion followed 
by periods of relatively quiet accretion and flux accumulation lead to oscillations in both 
magnetic flux and jet efficiency (Figure 6f,g). In this MAD accretion regime, first simulated 
in the non-relativistic context (Igumenshchev, Narayan & Abramowicz 2003; Igumenshchev 
2008), the magnetic flux on the black hole is as strong as possible and is dynamically impor- 
tant. Its time-average value, (emp) ~ 50, translates due to eq. (14) into the time-average jet 
efficiency, (n) ~ 140%. That (7) > 100% implies that the black hole releases more energy 
in the form of jets than the entire energy it receives from the accretion flow. Of course, 
the total energy is conserved, and the extra 4096 of energy comes from the black hole spin: 
Black hole rotation slows down woing to the action of jets. In fact, because of this, in the 
MAD state black holes slow down to essentially zero spin (Tchekhovskoy 2015). 

So far, we have discussed jets from geometrically thick disks with aspect ratio h/r ~ 
0.3—1. Radiatively efficient disks, such as those thought to power luminous quasars, are 
much thinner, with h/r ~ 0.01(L/0.1Lg). Unexpectedly, such thin disks led to the pro- 
duction of jets with efficiency n = 20—50% (Liska et al. 2019c), defying the theoretical 
expectations that geometrically thin disks are incapable of dragging large-scale poloidal 
magnetic flux inward and powering the jets (Lubow, Papaloizou & Pringle 1994). In these 
simulations, the thin disk was formed by rapid cooling of a thick disk, which is similar to 
the physics expected to take place during the hard-to-soft spectral state transition. Thus, 
these powerful jets might represent the transient jets seen in XRBs. They could also be 
related to radio-loud quasars that make up 1096 of all quasars. It is possible that radio 
quasar is a transient stage in the life of a quasar, and the radio quasars are undergoing a 
spectral state transition from low-luminosity AGNs to radio-loud quasars. The simulations 
show signs that these jets are indeed transient phenomena:Tthe magnetic flux strength on 
the black hole appears to undergo steady decline, indicating that the large-scale magnetic 
flux diffuses out of the black hole and out to larger radii through the disk, suggesting an 
impending shutoff of the jets. However, longer-duration simulations are needed to verify 
this hypothesis. 


5.2. Origin of Large-Scale Vertical Magnetic Flux 


An important question arises: Do black holes in nature receive enough large-scale magnetic 
flux to produce such powerful jets? Indications are that such strong magnetic fluxes and 
powerful jets are present in a wide range of systems ranging from AGNs (Zamaninasab 
et al. 2014; Ghisellini et al. 2014; Nemmen & Tchekhovskoy 2015) to tidal disruption events 
(TDEs; Tchekhovskoy et al. 2014) and gamma-ray bursts (Tchekhovskoy & Giannios 2015). 
What is the origin of this magnetic flux? Whereas in AGNs the ISM contains a substantial 
amount of large-scale poloidal magnetic flux, sufficient to flood the black hole if the accretion 
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Figure 6: Large-scale vertical magnetic field accumulates at the center and forms a dynami- 
cally important magnetic field that obstructs the accretion flow and leads to a magnetically 
arrested disk (Tchekhovskoy, Narayan & McKinney 2011). Color shows the logarithm of 
density (see color bar). (a) The horizontal and vertical slices through the initial condition, a 
radially extended (from rin = 15rg out to rout ~ 2x 10^r,) equilibrium hydrodynamic torus 
embedded with a weak (plasma 8 = 100) poloidal magnetic flux loop, shown with directed 
black lines. Even though the magnetic field is not important dynamically, the large size 
of the torus allows it to hold a large amount of poloidal magnetic flux. (b) Gas brings in 
poloidal magnetic flux, which accumulates around the black hole and squeezes the accretion 
disk vertically. (c,d) The black hole periodically ejects excess magnetic flux into the disk in 
the form of low-density magnetic flux eruptions. (e) Mass-energy accretion rate on the black 
hole versus time. (f) Whereas initially the magnetic flux on the black hole monotonically 
grows in time, at t ^ 6000r,, it becomes dynamically important and partially leaves the 
black hole, leading to a drop in dn. Soon thereafter, the system settles into a quasi steady 
state: The black hole overeats magnetic flux before shedding its excess into the accretion 
disk, after which the cycle repeats again. (g) Jet energy efficiency undergoes oscillations 
that mirror those of opu. The time-average value of the efficiency shown with the dashed 
line, 7 ~ 140%, exceeds 100%. This is the first demonstration of net energy extraction from 
an accreting black hole in a numerical simulation. 


disk can drag the flux toward the black hole (Narayan, Igumenshchev & Abramowicz 2003), 
the origin of the flux is less clear in TDEs when an unlucky star wanders too close to a 
supermassive black hole and gets spaghettified by its tidal forces. This is thought to result 
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in an accretion disk that feeds the black hole for months to years. In one remarkable TDE, 
the formation of the putative disk was accompanied by the launching of a powerful jet 
(Bloom et al. 2011; Burrows et al. 2011; Levan et al. 2011; Zauderer et al. 2011) that would 
require orders-of-magnitude-more large-scale flux on the black hole than that available in the 
disrupted star (Tchekhovskoy et al. 2014). One possibility is that the flux could be provided 
by a preexisting accretion disk (Kelley, Tchekhovskoy & Narayan 2014). However, could 
the large-scale magnetic flux be produced in situ, by the turbulent dynamo in the accretion 
disk? The problem of large-scale poloidal flux dynamo is long-standing and interdisciplinary: 
'The causes the of 11-year cycle of magnetic flux polarity flips on the Sun and the much 
longer dynamo cycle in Earth's core (Cowling 1981) are still unclear. Numerical simulations 
indicate that the formation of powerful jets requires preexisting large-scale vertical magnetic 
flux in the flow (Beckwith, Hawley & Krolik 2008; McKinney, Tchekhovskoy & Blandford 
2012). More recently, evidence points to the large-scale poloidal flux dynamo producing 
large-scale magnetic flux at least in some circumstances; in fact, in thick disks, h/r ~ 0.3, 
the generated magnetic flux becomes so strong that it leads to the formation of a MAD, 
even in the absence of any vertical magnetic flux to start with (Liska, Tchekhovskoy & 
Quataert 2020). However, thinner disks h/r ~ 0.02, so far do not show evidence of this 
process (Liska et al. 2019a). 


5.3. Jet Acceleration 


Many AGN jets reach relativistic velocities, with typical values of bulk Lorentz factor of 
^ ~ 10. Relativistic motion is not limited to AGNs, but is also observed in XRBs, with 
y ~ few, and gamma-ray bursts, with y 2 100. How do jets accelerate to such remarkably 
high velocities? They do so by converting their internal (magnetic, thermal) energy into 
bulk kinetic energy (Li, Chiueh & Begelman 1992; Beskin & Nokhrina 2006). Typically, 
such acceleration is accompanied by collimation. For instance, in the case of a well-studied 


jet in the M87 galaxy, the Lorentz factor is observed to increase as the jet half-opening angle 
—0.423-0.02 


decreases as a power law in radius, 0 xr , over more than five orders of magnitude 
in distance before levelling off at 0 ~ 0.01 rad ~ 0.6 degrees (Nakamura & Asada 2013), 
and the jet reaches Lorentz factors of y ~ 6 (Mertens et al. 2016). In idealized settings, 
the jet Lorentz factor increases inversely proportional to its half-opening angle, such that 
y ~ 67" (Tchekhovskoy, McKinney & Narayan 2008). The physics behind this acceleration 
law is simple: The jet accelerates such that its opening angle 0 does not exceed its beaming 
angle, » !, so the jet maintains transverse causal contact, so that the interior of the jet 
avoids running into the edge of the jet. However, this acceleration cannot last indefinitely: 
Eventually, the jet converts most of its internal energy into bulk motion kinetic energy, and 
the acceleration levels off at y ~ Ymax. Recently, Chatterjee et al. (2019) simulated the 
propagation of black hole jets over five orders of magnitude in distance and showed that the 
development of the single-power law jet shape extending for as many orders of magnitude 
occurs naturally if the jets are collimated by an extended accretion flow. Furthermore, the 
jet accelerates in a way similar to observations. 


5.4. Jet Emission 


It is broadly agreed upon that jet emission is due to a combination of synchrotron and inverse 
Compton processes. However, presently, there is no agreement on what accelerates the 
emitting electrons. In particular, we do not understand what accelerates the electrons that 
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produce the radio emission that we use as an indicator of jet activity. The problem is even 
more severe in the origin of high-energy gamma-ray (GeV and TeV) emission, where it is 
unclear whether the emission is coming from near the event horizon or from large distances. 
There are several candidate processes potentially responsible for accelerating the emitting 
electrons and, possibly, positrons. They can be accelerated in magnetospheric gaps near 
the black hole event horizon (e.g., Hirotani & Okamoto 1998; Broderick & Tchekhovskoy 
2015; Hirotani & Pu 2016; Ptitsyna & Neronov 2016; Levinson & Segev 2017; Chen, Yuan & 
Yang 2018; Parfrey, Philippov & Cerutti 2019), in 3D magnetic kink instabilities and shocks 
that can be triggered by jets running into the ISM (e.g., Bromberg & Tchekhovskoy 2016; 
Tchekhovskoy & Bromberg 2016; Barniol Duran, Tchekhovskoy & Giannios 2017), and in 
interface instabilities (e.g., Chatterjee et al. 2019). Heating due to magnetized turbulence 
in the exterior to the jet is also possible, so what we are seeing could be the emission from 
the sheath that surrounds the jet (and not from the jet itself Moscibrodzka & Falcke e.g., 
2013; Ressler et al. e.g., 2017). 


5.5. Tilted Disk Precession, Jets, Alignment, and Tearing 


A standard approach to modeling black hole accretion is to consider an accretion disk 
lying in the equatorial plane of the black hole. However, typically the disk midplane is 
tilted relative to that of the black hole, necessitating the consideration of tilted disks. Early 
simulations of tilted thick disks with h/r ~ 0.2 (Fragile & Anninos 2005; Fragile et al. 2007) 
confirmed the analytical expectations that general relativistic frame dragging by spinning 
black holes causes the tilted disks to undergo solid-body precession. However, what happens 
to their jets? Do they point along the black hole spin axis or disk rotation axis? The jets 
fly out along the direction of the disk rotation axis (McKinney, Tchekhovskoy & Blandford 
2013; Liska et al. 2018). However, if the magnetic field is dynamically important, the black 
hole manages to bring the disk and the jet into alignment at small radii, with the jet initially 
flying out along the black hole rotational axis before aligning with the rotational axis of the 
disk at large radii (McKinney, Tchekhovskoy & Blandford 2013). As the disk precesses, the 
large-scale jet precesses as well (Liska et al. 2018, 2019b), enabling the use of precessing 
jets as probes of strong-field gravity and general relativistic frame dragging. 

The response to the tilt of thinner disks is qualitatively different than it is to that of 
thick disks. In fact, when the disk thickness is smaller than the disk viscosity parameter, 
h/r « a, warps propagate in the disks viscously. Bardeen & Petterson (1975) predicted 
that this would lead to the inner parts of the accretion disk aligning with the black hole 
midplane and separating from the outer, misaligned part of the disk by a smooth warp (the 
Bardeen-Petterson effect or BP). However, at the time it was impossible to take into account 
the magnetized turbulence responsible for the accretion. And, though such alignment was 
seen in non-elativistic smoothed particle hydrodynamics (SPH) simulations (e.g., Nelson & 
Papaloizou 2000), until recently general relativistic magnetized simulations of tilted disks, 
as thin as h/r — 0.08, have shown no sign of the BP alignment (Zhuravlev et al. 2014). Liska 
et al. (2019c) simulated the thinnest disk, with h/r — 0.03, tilted by 10 degrees relative to 
a rapidly spinning black hole, and found that the inner regions of the disk, r S 5r, aligned 
with the black hole. This remarkable discovery is the first demonstration of BP alignment 
in a general relativistic numerical simulation of a magnetized accretion disk. Importantly, 
the duration of these simulations is long enough for the inner regions of the accretion disk 
to achieve inflow equilibrium and quasi-steady state. The BP alignment was also seen in 
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shorter nonrelativistic simulations that approximated the effects of black hole rotation via 
additional torques in MHD equations (Hawley & Krolik 2019) and did not include the 
effects of apsidal precession, which can qualitatively affect the alignment process (Nealon 
et al. 2016). 


Figure 7: Tilted thin accretion disks tear up into individual, independently precessing 
subdisks. Whereas at small radii the jets align with the inner subdisks, at large radii 
the jets align with the outer subdisk. As the disks precess, jets run into subdisks that 
get in the way and can not only expel them, as seen in the movie of the simulation 
(https://youtu.be/mbnG5.UTTdk), but also lead to energy dissipation that modifies the 
radial emission profile of the disk. Figure adapated from Liska et al. (2020). 


Interestingly, a disk tilted by 65 degrees got torn up into several individually precessing 
subdisks (Liska et al. 2020), as seen in Fig. 7. Note that tilted disks in SPH simulations also 
show tearing but with an important difference: Tilted disks get torn into a large number of 
thin rings (Nixon et al. 2012). This difference likely stems from the effects of magnetized 
turbulence and large-scale magnetic fields that hold the disk together differently than hy- 
drodynamic viscosity. The observational manifestations of torn disks are far reaching. The 
interactions between adjacent disks can lead to additional dissipation and emission at large 
radii. Jets and radiation from the inner disk interacting with outer subdisks can contribute 
additional heating. Not all gas flows between the adjacent subdisks, and an interesting frac- 
tion blown away owing to the interaction with the jets and winds. These factors can lead 
to modifications in the emission profile, potentially reducing the tension with the observed 
disk sizes (see Fig. 1). Disk tearing can also lead to a wide range of variability, for example, 
sub-disks precessing through our line of sight might act as absorbers, from time to time 
dimming the emission of the central source. 


5.6. Jet Interaction with Ambient Medium 


As the jets emerge from the black hole's sphere of influence, they run into the ISM. At 
this point, their behavior qualitatively changes. In the case of the M87 jet, the jet stops 
collimating (Section 5.3) and starts to decelerate (Mertens et al. 2016). Jet interaction 
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Figure 8: (a) Low-power AGN jets (blue-green) succumb to global magnetic instabilities, 
stall within their host galaxies, and inflate quasi-spherical cavities (yellow). (b) High-power 


jets maintain their stability, leave their host galaxies, and form strong backflows. Thus, 
magnetic instabilities can be the key to resolving a 40-year long puzzle on the cause of 
Fanaroff & Riley (1974) morphological dichotomy of AGNs (Tchekhovskoy & Bromberg 
2016). 


with the ISM is poorly understood. For instance, there is no agreement on the origin of 
the Fanaroff & Riley (1974) morphological dichotomy of AGN jets: Fanaroff-Riley type 
I (FRI) jets appear to develop instabilities early on and often disrupt inside the galaxy, 
whereas FRII jets appear well collimated and stably propagate to outside of the galaxy. 
'This can have important consequences for their parent galaxy: Whereas FRII jets leave 
the galaxy unscathed and deposit their energy outside the galaxy, FRI jets inject their 
energy into the galaxy and can significantly affect the star formation and dynamics of gas. 
'Thus, it is important to understand the stability properties of the jets. Whereas it is rather 
easy to reproduce the stable FRII jet morphology in numerical hydrodynamic simulations 
(Clarke, Norman & Burns 1986), reproducing the FRI morphology turned out to be much 
more difficult. Some of the possibilities include Kelvin-Helmholtz (KH) instabilities in the 
shear layers (e.g., Kaiser & Alexander 1997; Perucho, Martí & Hanasz 2005; Meliani & 
Keppens 2009; Perucho et al. 2010), and mass entrainment from stellar winds (Komissarov 
1994; Perucho et al. 2014; Wykes et al. 2015). Studies of magnetized jet stability were 
more promising; however, it turns out that magnetized jet stability sensitively depends 
on the degree of azimuthal winding of the magnetic field (e.g., Guan, Li & Li 2014), a 
free parameter whose value is poorly understood. This introduces an uncertainty into the 
factors that control jet stability. Tchekhovskoy & Bromberg (2016) carried out large-scale 
simulations of magnetized relativistic jets interacting with the ISM that for the first time 
were launched via the rotation at the base, as they are launched in nature. The advantage 
of this approach is that it organically determines the degree of azimuthal winding and 
therefore leads to the same stability properties of the jets as in nature. Figure 8 shows 
that a change by two orders of magnitude in jet power leads to a drastic change in jet 
morphology: More powerful jets are stable, whereas less powerful jets are unstable, and 
this is in qualitative agreement with the FRI/II dichotomy. Future work including magnetic 
fields and jet precession, which naturally emerges as a result of tilted accretion, will help to 
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refine the models of jet feedback so that they can be used as subgrid models in cosmological 
simulations. 


6. SUMMARY AND OUTLOOK 


Our understanding of the central engines of AGNs has greatly advanced over past several 
decades. We have a fiducial model for how accretion proceeds in the general relativistic 
spacetime of a supermassive black hole and gained an understanding of the central role 
played by magnetic fields in the process of jet formation and angular momentum trans- 
port. This general picture of accretion-powered emission is well supported by the existing 
observational constraints and we owe many of these insights (in part) to the application of 
state-of-the-art numerical simulations of GRMHD flows. This is particularly true for our 
understanding of jet production and jet-disk interactions. We are beginning to understand 
the effects of magnetic field geometry on the dynamics of accretion flows, the ability of disks 
to transport large-scale magnetic fields, and the disk-jet connection. The jets appear to be 
more resilient than previously thought and defy the standard expectations: They emerge 
in systems (a) without any large-scale vertical magnetic flux, (b) with a large disk tilt, and 
(c) with extremely small disk thickness. 

However, many challenges remain with our theoretical understanding of accretion disks. 
Longstanding theoretical inconsistencies remain unresolved, and detailed comparisons be- 
tween observations and theory yield significant discrepancies. What many viewed as promis- 
ing early agreement between theory and observations has faltered in the face of new ob- 
servational constraints: far-UV SEDs, size-scale constraints, and variability studies. For 
these reasons, we have chosen to focus much of this review on the limitations of our current 
understanding, with the hope of providing a framework to motivate future research. We 
have argued that global numerical simulations will be essential for solving some of the most 
important theoretical issues, including the following: 


e It is possible that gas fed to the accretion flows on large scales may be misaligned with 
the black hole spin. We must understand the warping of the accretion flow in such 
cases, including the possible impact of disk tearing, and its impact on disk emission 
and reprocessing. 

e Standard accretion disk models predict the radiation pressure-dominated central re- 
gions of accretion flows in AGNs to be violently unstable to thermal and inflow 
instabilities, yet there is no clear evidence of instability in the generic variability of 
these sources. Local MHD simulations cannot address inflow instability and provide 
conflicting results on the presence of thermal instability. 

e Local simulations and simple models suggest that opacities due to atomic transitions 
could strongly modify the structure of AGN accretion flows or even drive outflows. 
Such dynamics are already seen in simulations of massive stellar envelopes. 

e Magnetic pressure support seems to be a promising candidate for stabilizing accretion 
flows in the radiation pressure-dominated regime. Such support also leads to thicker 
accretion flows and both may change the vertical structure of accretion flows and 
increase the reprocessing of emission from the inner regions of the flow. 


AII the above questions require global simulations with large dynamic range. They may also 
require resolving relatively thin accretion flows, making them computationally expensive. 
Many of these questions also require increasingly complex treatments of radiation transfer 
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built on top of already sophisticated algorithms for evolving GRMHD. Progress will almost 
certainly require utilizing efficiently optimized algorithms that can harness the capabilities 
of the next generations of computing infrastructure (exascale and beyond). Fortunately, 
there has already been substantial recent progress in developing the numerical tools and 
algorithms needed for the future, and we are optimistic these questions can be addressed 
in the next decade. 


FUTURE ISSUES 


1. Although local simulations have advantages for studying accretion disks at high res- 
olution, questions about convergence and dependence on simulation domain sizes 
and magnetically powered outflows strongly motivate global numerical simulations 
of accretion flows. The large dynamic range required by such simulations requires 
the development of algorithms and codes that model the required physics while effi- 
ciently utilizing the latest advances in high-performance computing infrastructure. 

2. Multiscale simulations that bridge the scales of the galaxy and that of the black 
hole will be needed to self-consistently determine the mass supply near the black 
hole and the feedback of the black hole on the galaxy. 

3. Realistic treatments of radiation transfer are essential for modeling the radiation- 
dominated regions of accretion flows to resolve longstanding questions about stabil- 
ity and to study the effects of atomic opacities. The ability of scaling simulations 
to different black hole masses that is present in most previous nonradiative simula- 
tions is no longer possible when radiation transfer is included, requiring AGN- and 
XRB-specific simulations. 

4. Future simulations will need to address the evident discrepancies between observa- 
tions of AGNs and standard disk models. This includes size discrepancies, flattening 
of spectra in the UV and absence of edges, the origin of AGN continuum variability, 
and lack of generic evidence of instability. 
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